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We use series expansion techniques for analyzing properties of the phase transition between the 
Mott insulating and superfluid phase for bosons on the kagome lattice, and the multicritical point 
in the ground-state phase diagram for unit-filling is calculated. It is seen that of the clusters that 
contribute with non-zero weights to the ground state energy, many contain rings. The exponential 
decay coefficients of ground state correlations are also obtained within the Mott phase. For excited 
properties, quasiparticle dispersion and effective masses for particles and holes are computed. Fur- 
thermore at 8**^ order, the coherence-length critical exponent i/ is found to be comparably close to 
that of the 3D XY model. 
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Strong coupling expansion has proved to be one 
among many powerful techniques for providing quantita- 
tive descriptions of strongly correlated bosonic systems. 
Even at low order—, the disctinction between ID and 
higher dimensional critical systems was clearly revealed. 
This was confirmed in a DMRG study^ and a high order 
series expansion^. Prior to this, an infinite range mean- 
field analysis^ had proposed the various possible phases 
in such systems and a later local mean-field study^ pro- 
vided a quantitative description of the phase diagram. 
Furthermore, the suggestion that ultracold dilute gases 
of bosonic atoms could be modelled by such systems^ 
and the subsequent experimental realization of the tran- 
sition from the Mott insulator (MI) to the superfluid (SF) 
phased, renewed interest in the study of the field within 
the context of cold atoms. 

The current work analyzes such bosonic Hubbard 
models on the kagome lattice (see Fig. [T]) using the se- 
ries expansion technique^, which has been applied on the 
2D square lattice^ and the 2D triangular lattice^. The 
present study attempts to quantify the bosonic model's 
characteristic features on the kagome lattice, comple- 
menting earlier studies of similar systemsiSr— . De- 
spite the enormous bosonic Hilbert space, the similarity 




(a) (b) 



FIG. 1: (a) Kagome lattice in real space with basis vectors 
(bold) and unit cell (dotted), (b) Reciprocal lattice with sym- 
metry points labeled. 



transformations* can be used to treat such systems by 
creating effective Hamiltonians; in addition, there are no 
finite size effects since the thermodynamic limit is taken 
care of right at the very outset of the calculations^. 

The motivation to study such systems on frustrated 
lattices can be traced back to the suggestionii that the 
resonating- valence bond is the ground state of quantum 
antiferromagnets. Though there are some properties of 
the ground state of such systems - like the absence of long 
range magnetic order even down to low-temperatures 
- that are agreed upon, questions regarding the low- 
lying excitations or the gaps are not fully resolved. It 
is likely that such questions can be resolved through the 
setup of atomic quantum gases on similar optical lattices. 
An experimental setup for constructing and controlling 
optical kagome lattices using standing waves has been 
proposedi^. Such experiments could give us a clearer 
understanding about the corresponding spin-models via 
the relevant mapping between spins and bosons^. 

The Hamiltonian that describes such a model is the 
Bose-Hubbard Hamiltonian (BHH)^, 

(1) 

where the foj and bi are bosonic creation and annihilation 
operators, lii = bjbi is the number operator, the hopping- 
terms t are between nearest neighbors, and the system 
consists of a single species of soft-core bosons. The local 
energy term U contributes to a repulsive on-site interac- 
tion between bosons and fj, is the chemical potential. 

Before we head onto the results, clarifying the mean- 
ings of a few technical terms is in order. A topological 
graph is a graph which is completely specified by its ad- 
jacency matrix, which is a representation of the intercon- 
nectedness of its vertices; a topological graph may have 
many distinct embeddings on the lattice (corresponding 
to various orientations), each of which is called a clus- 
ter. It is only required to choose one representative from 
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FIG. 2: (Color online) Ground state energy of the Bose- 
Hubbard model in Equation (1) on the kagome lattice. 

these clusters (for a given topological graph) while the ge- 
ometric contributions of the clusters on the lattice can be 
encapsulated in an appropriate lattice constant. A uni- 
cyclic (bicyclic) graph contains one (two) cycles within 
the graph whereas a tree is acyclic. In the rest of the pa- 
per, we call unicyclic and bicyclic graphs as ring graphs. 

I. GROUND STATE PROPERTIES 

In this section, we present our results on a system- 
atic strong coupling expansion of ground state properties. 
Firstly, we calculate the ground state energy per site of 
the system in the thermodynamic limit, to 8*^ order in 
perturbation theory as 

^ = -8A2 + 24A3 - 40A4 + 
N 3 

28640 6 1706768 7 7644904 g 

where A = —t/U and N is the number of lattice sites. 
The expansion is obtained by starting the system at unit- 
filling and applying the hopping perturbation to the given 
order. 

From embedding considerations on the lattice, 152 
topological graphs contribute to the above expression for 



Graph c 


Starting contribution S{c, C) 


P2 


-8A^ 


Cs 


24A^ 


Ps 


-72A* 


^3.1 


677.33A'' 


Butterfly 


-2465.78A'^ 


P4 


-1777.78A'' 



TABLE I: Graphs c with starting contribution S{c, C) to 
weighted weights L{c,C)W{c) > L(C6,C)W{Ce) to the 
ground state energy in Equation (2). 



an 8*^ order calculation on the kagome lattice. And from 
these 152 graphs, only 19 graphs have non-zero weights 
W{c) associated with them; W{c) is defined as 

El^J2Lic,C)W{c) 

c 

where the summation is over all topological graphs c, 
L{c, C) is the lattice constant of the graph c on the lattice 
C and the W{c) is the weight (a corresponding polyno- 
mial) of the graph c. For bosonic hopping models, the 
lowest order of a contributing graph would be the short- 
est path on the graph to make a tour; this was noted 
previously2£, and the problem can be reformulated into 
the well-known Chinese Postman problem in graph the- 
ory. The validity of this optimization has been checked 
for ground state energy and effective Hamiltonian calcu- 
lations of other lattices as well. For the same order, we 
compare these with the ground state calculations of the 
triangular and square lattices in 2D: only 38/321 topo- 
logical graphs and 18/125 topological graphs contribute 
to these two lattices respectively, after the optimization. 

We further notice that of the 19 graphs, 12 of 
them have closed ring structures (unicyclic and bicyclic 
graphs); we isolate their contributions from the trees in 
Fig. [2l In the terminology of modern graph theory, the 
ring graphs are 3 tadpole graphs (Ta,!, 2, Te^i), 4 cycle 
graphs (C3, Cg, C7, C's), a buUgraph, a butterfly graph, 
a 3-bellgraph, a cricket graph and graph Xiq^t^^^- The 
trees are the 4 path graphs (P2, -P3, ^4, ^5) and 2 star 
graphs (iSs, 5*4). Thus, it is conceivable that new and in- 
teresting physics may be seen in the ground state prop- 
erties of bosonic or spin models on the kagome lattice 
by including more than hexagonal-interaction graphs in 
the proposed ring-exchange model^^i^. For example, 
the hexagon graph has, including the lattice constant, a 
maximal starting contribution of — 1576A^ to the ground 
state energy. Graphs which have higher contributions are 
listed in Table U 

In this ground state, we further calculate the corre- 
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0.006 
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0.0049 


0.018 
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0.0095 


0.024 
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0.015 
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0.034 
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0.048 


1.53 


0.058 
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1.35 
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1.17 
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TABLE II: Exponential decay coefficients Oft/jy for the cor- 
relator in Equation (3) with their corresponding errors, for 
varying perturbation strengths. 
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FIG. 3: (Color online) (a) Quasiparticle dispersion for particle and hole along a given symmetry line in the first Brillouin zone 
for t/U — 0.050. (b) Effective mass tensor component mn for particle and hole excitations near the bottom of the band. The 
point F is defined to be approximately one-tenth of the distance along F — M. The symbols indicate values obtained from 
numerical differentiation of a 15*'^ order Chebyshev-approximated dispersion tk as in Equation (4) . 



lator 



(3) 



and determine the exponential decay coefficients i.e. 

L ) . We see from Table |TI] that as we 
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restore the U{1) symmetry and enter deeper into the MI 
phase, the boson particle-number spatial correlations are 
captured by the current method with decreasing errors. 



II. EXCITED STATE PROPERTIES 

Single particle (hole) excitations are created within 
the insulating phase and their low-lying dispersion curves 
arc obtained using the similarity transformation. The ex- 
citations on the lattice A(r), obtained by the embedding 
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FIG. 4: (Color online) (a) Extrapolation of the zeros of the 
truncated gap-polynomial in Equation. (5) as a function of 
inverse order; the errors in the plotted zeros of the finite poly- 
nomials are ignored. 



of the effective Hamiltonians on the lattice, are diagonal- 
ized by plane-waves^ 



e(k) 



'A(r), 



and the spectra can be obtained for a given value of the 
perturbation strength. For t/U ~ 0.050, the particle and 
hole spectra are shown in Fig. |3(a) The smallest value of 



the particle-hole gap can be clearly seen at the k = (0,0) 
point. We also notice a nearly flat band in the M ~ X 
direction in k-space. 

Once we obtain the dispersion within this formal- 
ism, a straightforward calculation can give us the effec- 
tive mass tensor components rriij for particles and holes 
as 



92e(k) 
dkidkj 



The mass tensor component mn, with ki = k^, for par- 
ticle and hole excitations are shown in |3(b)| for the same 
perturbation strength along a given direction in the first 
Brillouin zone (FBZ) near the bottom of the band. The 
solid lines indicate the tensor component calculated di- 
rectly from the excitations A(r). The symbols in the 
figure arc obtained as follows: a 15"^ order Chebyshev 
approximation of eyr^ is obtained as 
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(4) 



where T„(k) is a Chebyshev polynomial of degree n, N ■ 
15 and the coefficient Cj is given by 
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The approximation in Equation (4) is numerically dif- 
ferentiated twice to give the inverse-mass tensor compo- 
nents. The variables have been suitably rescaled to fit 
the range [—1, 1] where the above approximation holds. 
Aawy from the bottom of the band, we observe recurring 
divergences in the mass components within the FBZ. It 
is unclear how these can be renormalized into physically 
valid results within the current formalism. 

Finally we arrive at an expression for the particle- 
hole gap (5(k = 0) in the insulating phase by pertur- 
batively calculating the particle and hole contours of the 
ground state phase diagram and taking the difference be- 
tween the two. This gives 



(5(k = 0) = 1 + 12A - 22A^ + 75A^ 
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(5) 



Closing of the gap will indicate the complete dis- 
appearance of the insulating phase. To this end, two 
approaches are undertaken to analyze the above series. 
First, we truncate the above polynomial to order n and 
solve for the zeros of the truncated polynomial and iden- 
tify the relevant zero as t", with U = 1. This is then ex- 
trapolated to an infinite order polynomial, which should 
directly give an estimate of the critical point. This gives 
the critical point as tc = 0.0617 ± 0.0001 as seen in Fig. 
II 

Secondly, since the gap obeys a power law of the 
forra^ 

s{k = o) = A{u-ty', 

where v is the coherence-length critical exponent and z is 
the dynamical critical exponent, we can perform a Fade 
approximant to the series in Equation (5) after taking its 
logarithm and differentiating it, to give 
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A{t) 



tr - t 



Considering the diagonal approximants ([2/2], [3/3], 
[4/4]) we find U, = 0.0615 ± 0.0001 and ly = 0.651 ±0.006 
after setting 2 — 1 at the multicritical point^. The 
critical point calculated by the direct extrapolation is 
seen to be a reasonably good estimate for an 8"^ order 



calculation; the critical exponent is to be compared to 
the corresponding value of the 3D XY model, which has 
1^ = 0.67155 ±0.00027-2^. It is interesting to note that for 
unit-filling on the kagome lattice, performing a mean- 
field calculation^! gives tc ~ 0.04289. For comparison, 
we compute the superfluid stiffness p using the ALPS 
software^S^ to perform a Quantum Monte Carlo simu- 
lation implemented by the worm algorithm for 3 different 
lattice sizes (number of kagome unit cells: 3x3, 6x6, 9x9), 
with a maximum of 2 bosons per site. Using the scaling 
ansat^ 



P 



we get an estimate tc ~ 0.061. Here, d is the dimension- 
ality, L is the linear system size, p is the scaling function 
and a is a non-universal metric factor. As before, we set 
z = 1 at the multicritical point. 

To conclude, we have seen that the number of con- 
tributing graphs to ground state energy calculations on 
the kagome lattice using the series expansion is only a 
small fraction of the number calculated from geometric 
considerations of embedding, which is a feature of hop- 
ping models. The optimization for such models, resulting 
from the use of the Chinese Postman algorithm to elim- 
inate graphs, has been checked. Furthermore, it is pos- 
sible that the inclusion of these ring graphs into the in- 
teraction terms of the previously proposed ring-exchange 
Hamiltonian might give further insight into the ground 
state properties of similar models defined on the kagome 
lattice. We also present an accurate determination of the 
critical point and critical exponent for the model along 
with other excited state properties, albeit the perturba- 
tive order was lower than previously studied 2D lattices; 
yet we believe that this is a reasonably precise estimate 
that can be used as a reference for future work. 
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